
a <- 3.06; b <- 2.56

n <-20; y <- 12

pred_p_sim <- rbeta(1, a+y, b+n -y)

pred_y_sim <- rbinom(1, n, pred_p_sim)


s <- 1000

prep_p_sim <- rbeta(s, a+y, b+n-y)
pred_y_sim <- rbinom(s, n, prep_p_sim)


hist(pred_y_sim, breaks = seq(-0.5, n + 0.5, 1),
     main = "Predicted Y from Beta-Binomial Simulation", xlab = "Predicted Y", ylab = "Frequency")

